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ABSTRACT 

New multi-roll coronagraphic images of the HD 181327 debris disk obtained using the Space Telescope Imaging 
Spectrograph on board the Hubble Space Telescope reveal the debris ring in its entirety at high signal-to-noise ratio 
and unprecedented spatial resolution. We present and apply a new multi-roll image processing routine to identify 
and further remove quasi-static point-spread function- subtraction residuals and quantify systematic uncertainties. 
We also use a new iterative image deprojection technique to constrain the true disk geometry and aggressively 
remove any surface brightness asymmetries that can be explained without invoking dust density enhancements/ 
deficits. The measured empirical scattering phase function for the disk is more forward scattering than previously 
thought and is not well-fit by a Henyey-Greenstein function. The empirical scattering phase function varies with 
stellocentric distance, consistent with the expected radiation pressured-induced size segregation exterior to the belt. 
Within the belt, the empirical scattering phase function contradicts unperturbed debris ring models, suggesting the 
presence of an unseen planet. The radial profile of the flux density is degenerate with a radially varying scattering 
phase function; therefore estimates of the ring’s true width and edge slope may be highly uncertain. We detect large 
scale asymmetries in the disk, consistent with either the recent catastrophic disruption of a body with mass >1% 
the mass of Pluto, or disk warping due to strong interactions with the interstellar medium. 
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1. INTRODUCTION 

Images of spatially resolved debris disks show a range of 

spectacular asymmetries including eccentric rings (e.g., Kalas 
et al. 2005; Schneider et al. 2009; Krist et al. 2012; Boley et al. 

2012), warps and sub- structures (e.g., Golimowski et al. 2006; 
Krist et al. 2005), and various other morphologies (e.g., Hines 
et al. 2007; Kalas et al. 2007). Interpreting surface brightness 
asymmetries as dust density variations and explaining them via 
dynamic processes has been popular since coronagraphic im- 
ages of the /3 Pictoris disk revealed scattered light asymme- 
tries (Kalas & Jewitt 1995). Recently, ALMA has begun to ex- 
plore such patterns at unprecedented resolution at submillimeter 
wavelengths (e.g., Boley et al. 2012). 

Many models have shown that exoplanets can create asym- 
metric dust distributions in debris disks via gravitational pertur- 
bations, potentially revealing the presence of otherwise unde- 
tectable planets. Disk asymmetries created by planets could be 
the only way to detect true Neptune analogs orbiting nearby stars 
on reasonable timescales. Direct images of exoplanet candidates 
associated with debris disks have begun to demonstrate both the 
potential and complexities of this concept for locating new plan- 
ets and constraining their properties (e.g., Quillen 2006; Chiang 
et al. 2009; Lagrange et al. 2010; Kalas et al. 2013). 

Many other dynamical processes can also produce dust 
density asymmetries in debris disks. Disks can interact with 
the interstellar medium (ISM; Artymowicz & Clampin 1997; 


Debes et al. 2009; Maness et al. 2009; Marzari & Thebault 
2011). In the solar system, collisions of asteroids can produce 
detectable trails of debris (Jewitt et al. 2010). In debris disks, 
recent collisions can potentially yield detectable arcs of debris 
(Grigorieva et al. 2007; Krai et al. 2013; Jackson et al. 2014). 

It remains crucial to pursue explanations for debris disk asym- 
metries other than density enhancements /deficits. The observed 
scattered starlight from a disk is a complex combination of the 
disk’s geometry, illumination, and stellocentric grain-size seg- 
regation, the dust grains’ size-dependent scattering efficiency 
and scattering phase function (SPF), and line-of-sight projec- 
tion effects in the case of inclined disks (e.g., ansal brightening 
for disks with large scale heights). All of these effects are de- 
generate to some degree and, given our current understanding 
of dust grain scattering properties, are exceedingly difficult to 
disentangle. 

In light of this, we interpret new multi-roll Space Telescope 
Imaging Spectrograph (STIS) coronagraphic images of the 
HD 181327 debris disk by aggressively pursuing explanations 
for observed asymmetries other than density enhancements. We 
exploit a symmetry in the SPF to search for deviations from a 
smooth disk. 

HD 181327 is an F5/6, ~12 Myr old main sequence member 
of the ft Pic moving group, located at a distance of 51.8 pc 
(Schneider et al. 2006). HD 181327 has a strong thermal IR ex- 
cess (L ir /L* = 0.25%) attributed to re-radiating circumstellar 
dust. This debris disk was first detected with IRAS (Mannings 
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Table 1 

Observations 
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Target 

Orientation 

(°) 

Exposures 

Total Exposure Time 
(s) 

Wedge 

Date 

HD 181327 

222.738 

8 

175.2 

0.6A 

2011 May 20 


222.738 

8 

175.2 

0.6A 

2011 May 20 


222.738 

4 

87.6 

0.6A 

2011 May 20 


222.738 

4 

1708.0 

1.0A 

2011 May 20 


242.738 

8 

175.2 

0.6A 

2011 May 20 


242.738 

8 

175.2 

0.6A 

2011 May 20 


242.738 

4 

87.6 

0.6A 

2011 May 20 


242.738 

4 

1708.0 

1.0A 

2011 May 20 

HD 180134 

243.219 

8 

95.2 

0.6A 

2011 May 20 


243.219 

8 

95.2 

0.6A 

2011 May 20 


243.219 

8 

95.2 

0.6A 

2011 May 20 


243.219 

8 

1656.0 

1.0A 

2011 May 20 

HD 181327 

262.738 

8 

175.2 

0.6A 

2011 May 20 


262.738 

8 

175.2 

0.6A 

2011 May 20 


262.738 

4 

87.6 

0.6A 

2011 May 20 


262.738 

4 

1708.0 

1.0A 

2011 May 20 


293.056 

8 

175.2 

0.6A 

2011 July 10 


293.056 

8 

175.2 

0.6A 

2011 July 10 


293.056 

4 

87.6 

0.6A 

2011 July 10 


293.056 

4 

1644.0 

1.0A 

2011 July 10 


313.556 

8 

175.2 

0.6A 

2011 July 10 


313.556 

8 

175.2 

0.6A 

2011 July 10 


313.556 

4 

87.6 

0.6A 

2011 July 10 


313.556 

4 

1708.0 

1.0A 

2011 July 10 

HD 180134 

314.232 

8 

95.2 

0.6A 

2011 July 10 


314.232 

8 

95.2 

0.6A 

2011 July 10 


314.232 

8 

95.2 

0.6A 

2011 July 10 


314.232 

8 

1656.0 

1.0A 

2011 July 10 

HD 181327 

334.056 

8 

175.2 

0.6A 

2011 July 10 


334.056 

8 

175.2 

0.6A 

2011 July 10 


334.056 

4 

87.6 

0.6A 

2011 July 10 


334.056 

4 

1708.0 

1.0A 

2011 July 10 


& Barlow 1998) and subsequently resolved at 0.6 /xm with 
the Advanced Camera for Surveys, 1 . 1 /xm with the Near In- 
frared Camera and Multi-Object Spectrometer on Hubble Space 
Telescope ( HST ; Schneider et al. 2006), 18.3 /xm with Gemini 
South T-ReCS (Chen et al. 2008), and 3.2 mm with the Aus- 
tralian Telescope Compact Array (Lebreton et al. 2012). All 
resolved observations are consistent with a ring of dust with 
radius ^90 AU and a cleared interior. Most resolved images 
suggest clumpy asymmetries, but also have low signal-to-noise 
ratios (S/Ns), so we refrain from incorporating those previous 
observations into our analysis. 

In Section 2, we briefly describe our new STIS observations 
of HD 181327 and present a new point-spread function (PSF) 
residual removal routine. For a more detailed description of 
the observations, see Schneider (2014). In Section 3, we detail 
our image deprojection techniques and ascertain the minimally 
asymmetric face-on optical depth. In Section 4, we interpret the 
observed SPF and discuss several explanations for the observed 
asymmetries. 


2. OBSERVATIONS AND DATA REDUCTION 

As part of the HST GO 12228 (PI: G. Schneider) observa- 
tion program, we observed HD 181327 in scattered light using 
the STIS coronagraph at six roll angles, each at two wedge- 
occulter positions (WedgeA0.6 and WedgeAl.0, which have 


occulting half- widths of O'.' 3 and (X'5, respectively 7 ). We ob- 
served HD 181327 over three orbits on 2011 May 20 and three 
orbits on 2011 July 10, and the (B — V color-matched) PSF 
reference star HD 180134 at each wedge position and a single 
roll angle interleaved with the target observations. Table 1 sum- 
marizes the observations. The STIS 50 CCD channel, equipped 
with the coronagraphic wedges used for our observations, has 
an image scale at the detector focal plane of 0'. / 05077 pixel -1 . 
Filters cannot be used with the coronagraphs, so the images are 
obtained with the full spectral response of the detector, i.e., a 
central bandpass of 0.5752 /xm and a FWHM of 0.433 /xm. 
Multiple exposures in each observational configuration (a given 
roll angle and wedge position) are median-combined. 

All astrometric measurements of the HD 181327 debris ring 
are referenced to the position of the coronagraphically occulted 
star, as measured individually from 12 independent images us- 
ing the diffraction spikes in each image to locate the star. In 
each of the 12 images, the uncertainty in the target position is 
approximately 0.3 pixels (0.8 AU) in the Science Aperture In- 
strument (detector) Frame (SIAF) v- andy-directions. In the 12- 
image combination, the uncertainty is reduced to ±0.087 pixels, 
or ±4.4 mas (0.23 AU); this is comparable to the HST point- 
ing stability (rms 2-guide star fine-lock jitter). The mean stel- 
lar position is used to anchor the inter- and intra- visit stellar 


7 See the HST STIS instrument handbook for a full description. 
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Number of observations 


Figure 1. Left: central 200 x 200 pixels of the 12-image median-combined image for HD 181327 (pixel scale = 50.77 mas = 2.63 AU). An artificial occulting 
spot with a radius of 18 pixels has been applied for illustrative purposes. Stellar position is marked with a star. Right: number of images used per pixel for median 
combination. 

(A color version of this figure is available in the online journal.) 


location following target acquisition slews and to reference the 
WedgeA0.6 pointings to the WedgeAl.O pointings. 

All visit-level PSF subtractions are done in the SIAF, treat- 
ing target and PSF-template positions and brightness as free pa- 
rameters while iteratively minimizing PSF-subtraction residuals 
following the procedure described by Schneider et al. (2009). 
Each image is then rotated about the mean stellar position to a 
common north-up orientation. We then manually build a mask 
specific to each observation, flagging those pixels that are ob- 
scured by the occulting wedge, corrupted by diffraction spikes, 
saturated or adversely affected by wedge-edge artifacts, or are 
beyond the field of view sub-array read out. Details of this 
process and considerations for optimization are discussed in 
detail in Schneider (2014). Finally, we median-combine the 
12 astrometrically co-registered images in the common celes- 
tial frame. The left panel of Figure 1 shows the inner 200 x 
200 pixels (pixel scale = 50.77 mas, or 2.63 AU assuming a 
distance of 51.8 pc) of our 12-image median. Our multi-roll ob- 
servation method reduces the mean inner working angle of the 
STIS coronagraph, reduces the influence of PSF artifacts, and 
improves S/N. 

The right panel of Figure 1 shows a map of the number 
of observations per pixel, n v i x . We mask off optical artifacts 
(occulting mask shadows, telescope diffraction spikes) that are 
rotationally invariant in the frame of the detector. Thus, with 
multi-image masking, each output pixel in the final image is 
derived from different numbers of input images, so not all pixels 
in the final image are identically exposed. The radial “spokes” 
in this map are primarily a result of masking off the telescope’s 
secondary mirror support structures at the six roll angles. For 
the WedgeA0.6 observations, the images become photon starved 
at a radius of ~2"6, such that the S/N ~ 1. To avoid adding 
unnecessary noise to the outer disk, we masked the WedgeA0.6 
images beyond 2" 6. As a result, the n v j x map exhibits a circular 
disk of radius 2"6, interior to which n pix is larger. 

Our multi-roll observation technique reduces both the impact 
of temporal instabilities in the PSF structures (“breathing”) and 
static PSF residuals that co-rotate with the instrument/ telescope 
optics. To estimate the remaining impact of PSF artifacts, we 
produced an 11 -image median, subtracted it from the 12-image 
median, and divided by the 12-image median to produce a 



< -5 -2.5 0 2.5 >5 

Difference (%) 


Figure 2. Relative difference between a representative 1 1 -image median and the 
12-image median (central 200 x 200 pixels). PSF residuals in a single image 
can bias the flux of the 12-image median by ~5% over scales ~10 pixels or 
larger. 

(A color version of this figure is available in the online journal.) 


fractional residual map. If the 12-image median is robust to 
these PSF residuals, then leaving out any one image should not 
greatly impact the final image. Figure 2 shows one such frac- 
tional residual map. The fractional residuals, smoothed using a 
3x3 pixel median boxcar and displayed on a saturated scale for 
illustrative purposes, show correlated biases in the median by 
as much as ~5% over scales ~10 pixels. Additionally, Figure 2 
shows that leaving out this particular image affects the NE-SW 
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<-10 -5 0 5 >10 


Relative flux (%) 

Figure 3. Static PSF residuals common to all images in the SIAF frame, 
smoothed by a 3 x 3 boxcar median (central 200 x 200 pixels). We remove 
these residuals from each individual image to create an improved multi-roll 
median image. 

(A color version of this figure is available in the online journal.) 

asymmetry at this level by enhancing the NE flux and reducing 
the SW flux. 

Although we cannot remove time-dependent PSF residuals 
with this method, we can better mitigate static PSF residuals. 
We did this using the following “multi-roll residual removal 
routine” (MRRR, pronounced “myrrh”): 

1. form the 12-image median, oriented north-up; 

2. subtract the 12-image median from each individual masked 
image, oriented north-up, to form 12 disk-less residual 
images; 


3 . de-rotate each disk-less residual image to the SIAF-oriented 
frame; 

4. take the median of the 12 disk-less residual images to 
produce a map of static residuals; 

5. smooth the static residual map using a 3 x 3 median boxcar 
to remove pixel-to-pixel noise, retaining only larger scale 
correlated residuals; 

6. subtract the smoothed, correlated residual map from each of 
the 12 images in the SIAF frame, rotate each to the north-up 
frame, and form a new residual-removed 12 image median. 

Figure 3 shows the smoothed, correlated residual map for the 
inner 200 x 200 pixels of our observations of the HD 181327 
disk (displayed on a saturated scale for illustration). Figure 4 
shows the final residual-removed 12-image median, and the 
fractional difference between the original 12-image median and 
the residual-removed 12-image median. In this case, MRRR 
dims the ansae by ~5% and brightens the SW side of the ring 
by ~5%. Not surprisingly, the regions of the disk impacted 
most significantly are the regions with small values of n p [ x (see 
Figure 1). We use the MRRR-corrected image for all subsequent 
analysis. 

One could in principle run MRRR several times over, iter- 
atively removing the correlated static residuals. A preliminary 
implementation of an iterative MRRR appears to converge in 
only a few iterations, with higher order iterations predominantly 
brightening the SW side of the disk by an additional ~2%. How- 
ever, MRRR requires the correlated PSF residual amplitudes to 
be larger than the Poisson noise remaining in the disk-less im- 
ages. It is unclear to us at what point this assumption breaks 
down, so for this work we conservatively perform only a single 
iteration. 

Modeling our MRRR-corrected image of the HD 181327 
debris disk requires an understanding of the uncertainty in the 
PSF-subtracted flux in each pixel. On average, the standard 
deviation of the flux measurements within each pixel greatly 
exceeds what is expected for Poisson noise alone, so the 
observations are dominated by systematic noise, as is typical 
for HST PSF template- subtracted imaging. In principle, one 
could use the standard deviation in each pixel as an estimate of 
each pixel’s uncertainty. However, with a median n p [ x = 6 near 
the ring’s peak flux, a large number of pixels will, by chance, 
have a standard deviation that is much smaller than the true 



0.00 1.48 2.96 4.45 5.93 -10.0 -3.3 3.3 10.0 

Flux density (counts second -1 pixel -1 ) Percent Difference (%) 

Figure 4. Left: final multi-roll image of HD 181327 after processing it with one iteration of MRRR. Right: fractional updates to the image produced by MRRR. 

(A color version of this figure is available in the online journal.) 
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systematic uncertainty. These pixels will dominate the error- 
budget of any model fitting routine. One could arbitrarily mask 
off pixels whose S/N exceeds some threshold, but that would 
remove useful data from the fit. 

As an alternative, we assume the dominant systematic noise is 
only a function of stellocentric angular distance. Because such 
noise is not dependent on disk brightness, we separately subtract 
each of our original 12 images from the MRRR-corrected image, 
then calculate the standard deviation of this set of differences 
as a function of stellocentric distance — this is equivalent to the 
standard deviation of the unsmoothed residuals. We then divide 
th e azimut hally symmetric standard deviation by the greater of 

y/n p ix — 1 and unity to account for the number of observations 
per pixel, n v j x , while avoiding a systematic underestimate of 
the uncertainty given small values of n p i x . Note that we base 

our uncertainty estimates on the first iteration residuals, which 
contain the PSF residuals removed by MRRR; our conservative 
uncertainty estimate budgets for the changes made by MRRR. 

3. RESULTS AND ANALYSIS 

The MRRR-corrected image of HD 181327, shown in 
Figure 4, exhibits several asymmetries that are immediately 
noticeable. The NE side of the disk exhibits a peak in the sur- 
face brightness, approximately 30% brighter than the SW side 
of the disk. There is also a NW-SE asymmetry, with the NW 
side of the disk approximately 10% brighter than the SE side. 
Additionally, the disk appears more radially extended toward 
the N than toward the S. 

It is not immediately obvious whether these asymmetries 
are due to local dust density enhancements or a combination 
of geometric and scattering effects. A belt with significant 
density enhancements may suggest the presence of a nearby 
perturbing planet — or that a massive collision recently occurred. 
Simulations of planets perturbing debris disks have shown that 
planets are capable of producing belts with sharp inner edges 
(e.g., Quillen 2006; Chiang et al. 2009; Rodigas et al. 2014) as 
well as azimuthal asymmetries in the dust density (e.g., Wyatt 
et al. 1999; Kuchner & Holman 2003; Ertel et al. 2012). 

The scattered-light asymmetries we observe in the 
HD 181327 disk are primarily located along the minor and ma- 
jor axes of the projected belt, suggesting that they may be due, 
at least in part, to geometric and scattering effects. Here we seek 
to determine whether geometric and scattering effects alone can 
explain these apparent morphological asymmetries, or if an ac- 
tual dust density enhancement is necessary. To do so, we assume 
that the disk is infinitely thin and flat (an assumption we address 
in Section 4.2), deproject the disk to a “face-on” viewing ge- 
ometry, and remove any asymmetries that can be explained by 
geometric and scattering effects. Any significant asymmetries 
that remain must therefore be density enhancements or deficits, 
or are signs that the disk is not flat. 

To deproject the observations and examine the face-on optical 
depth of the disk we employed the following process: 

1. fit the observed (projected) HD 181327 dust belt with an 
ellipse; 

2. deproject the ellipse to obtain the true orbital ellipse and 
disk geometry; 

3. using the true orbital ellipse values in the projected image 
plane, correct for the 1/r 2 illumination factor; 

4. determine the best-fit SPF and divide it out of the image; 

5. deproject the image to produce a final face-on optical depth 
image; 


6. examine the optical depth for remaining density asymme- 
tries. 

Below we discuss this process step by step. In practice, steps 
1-3 were incorporated into their own iterative subroutine which 
we describe below. We will repeatedly refer to Figure 6, which 
illustrates each step of this process. 

3.1. Projected Ellipse Fitting 

There are many ways to fit an ellipse to the HD 181327 dust 
ring shown in Figure 4. First, we must choose a metric that 
defines the ellipse, e.g., the location of the belt’s peak flux, an 
isophote, the belt’s inner edge, etc. Ideally we would choose a 
metric that reflects the underlying density distribution, e.g., the 
peak density of the belt. However, the radial location of the peak 
surface brightness does not correspond to the radial location 
of the peak density, because the differential 1/r * 2 3 4 5 illumination 
factor shifts the apparent peak location (an especially important 
factor for eccentric debris rings). We must correct for the 1/r 2 
illumination factor to fit the orientation of the disk, but we 
must know the orientation of the disk to correct for the 1/r 2 
illumination factor. 

In light of this catch-22, we developed a method to fit 
the illumination-corrected image by iterating steps 1-3 of our 
disk deprojection procedure. First, we guessed the orientation 
of the disk (e.g., circular and face-on). We then generated a 
map of 1/r 2 and used it to remove the stellar illumination. 
We then recorded the coordinates of the illumination-corrected 
projected belt maximum, fit these coordinates with an ellipse, 
and deprojected the ellipse to obtain a new disk orientation. 
We iterated this process, using updated 1/r 2 maps each time. 
It is unclear whether this procedure should converge in all 
cases, but in the case of HD 181327, we found that this process 
converged in ~3 iterations regardless of our initial guess of the 
disk orientation. 

We used polar coordinates to measure the projected radius 
of the peak illumination-corrected surface brightness, which we 
refer to loosely as the surface density. We calculated on-sky p 
and 0 values for each pixel’s center relative to the star location. 
We divided the surface density image into 172 wedges with 
angular size of 2°. 1 centered on the star, chosen such that the arc 
length of a wedge is ~1 pixel at the location of the projected 
belt’s semi-minor axis. 

For each wedge, we measured the sub-pixel radial location of 
the ring maximum by fitting a Gaussian to the radial profile near 
the peak surface density. To find the ideal number of pixels to fit, 
we made separate Gaussian fits to the nearest 7-17 pixels and 
chose the fit with the most certain peak location. This method 
produces a set of p pea k, the best-fit radial surface density peak, 
at each angular location, 0 pea k- 

We note that the uncertainty of the peak location of each 
wedge ultimately controls the uncertainty in the belt’s semi- 
major axis, eccentricity, and orientation. Thus, it is important 
to estimate the uncertainty in the peak robustly. To do this, we 
used a Monte Carlo procedure. Each time we fit the data with a 
Gaussian, we performed 100 Monte Carlo trial fits. For each of 
these 100 Monte Carlo trials, we added a random flux to each 
pixel, drawn from a normal distribution with standard deviation 
equal to the pixel’s flux uncertainty, and refit with a Gaussian. 
The uncertainty in the peak location was then set equal to the 
standard deviation of the Gaussian peak locations from the 100 
Monte Carlo trials. 

We then ran through a fine grid of ellipse parameters to search 
for the best fit to the (p pea k, 0 pea k) coordinates. For each set of 
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Table 2 

Best Fit Ellipses to the Projected HD 181327 Debris Belt 


Method 

a' 

(pixels)/(AU) 

e' 

PA' 

(°) 

Ax' 

(pixels)/(AU) 

Ay' 

(pixels) /(AU) 

Maximum 
Inner edge 

34.4^ 0 - 4 4 /90.5+ 1 1 1 i 

31.3^o 4 4 /82.3+ 1 i 1 i 

0 . 48 $% 

0 - 50 $% 

101.2^ 

98 - 5+40 

- 0 -34+ 0 6 3 3 5 8 /-0- 8 9- 0 i 9 o 2 o 

— 0.28+° 0 3 3 4 2 /-0.74+° 0 f 4 

0 52 +0 - 30 /I T7+ 0 - 79 
-0.30/ —0.79 

—0 1 1 +0 - 26 /— 0 29 +0 - 68 
U - 11 -0.24/ U ' Zy -0.63 


Table 3 

Deprojected Ellipse Parameters for the HD 181327 Debris Belt 


Method 

a 

(AU) 

e 

CO 

(°) 

i 

(°) 

Q 

(°) 

Maximum 
Inner edge 

90 - 5 - 1 i 1 i 

82.3+\' 1 i 

0 O2 +0 01 
u - uz -o.oi 

Q A 1+0.01 
u * ui -0.01 

-70+1 

16+\ 6 9 4 6 

28-5+2.0 

30-3+2X1 

ll-2! 4 4 6 6 

8-5+4X, 


ellipse parameters we calculated a x 2 value by comparing each 
(Ppeak ? 0 P eak) point with the nearest point to the ellipse, i.e., we 
minimized the perpendicular distance to each point weighted 
by the uncertainty. Table 2 lists the parameters of the best-fit 
ellipse to the illumination-corrected belt maximum. Here the 
primed quantities refer to the on-sky, projected ellipse. We fit 
the semi-major axis, a ', eccentricity, e ' , position angle, PA', 
and the location of the ellipse center (A*', Ay') relative to the 
star. To determine the la uncertainties, we normalized the x 2 
values such that the minimum x 2 was equal to unity, then took 
the projection of the parameters for all models with normalized 
X 2 < 2. We consider all projected ellipse fits with normalized 
X 2 < 2 to be acceptable fits, and analyzed all of these fits in the 
following sections. 

To check whether our results depended strongly on the ellipse 
metric we chose, we repeated this procedure by fitting the inner 
edge of the belt. To fit the inner edge, we first calculated the radial 
derivative of the illumination-corrected HD 181327 image in the 
true disk plane, for which the maximum marks the belt’s inner 
edge. We then fit the maximum of the radial derivative using 
the techniques described above. Table 2 lists the results of fits 
to the inner edge of the belt. With exception to the values of 
a' , which should not agree, and Ay', the fits agree to within la, 
suggesting that the disk is thin and flat near the brightest part 
of the observed ring. The small discrepancy in the Ay' values 
may be a result of the blurring of radial features along the disk 
minor axis due to a small non-zero disk scale height, which we 
address in Section 4.2. 

3.2. Ellipse Deprojection 

To deproject the ellipse fits, i.e., obtain the true geometry 
and orientation of the disk, we used the Kowalsky method, as 
described by Smart (1930). This analytic method transforms 
the parameters describing an apparent ellipse with a center 
offset from the star ( a ', e' , PA', Ax', Ay') into a set of unique, 
deprojected ellipse parameters (a, e, co, i, £2) describing both 
the geometry and orientation of the true ellipse. Here a is the 
true semi-major axis, e is the true eccentricity, co is the argument 
of pericenter, which defines the axis of inclination in the plane 
of the disk, i is the inclination, and £2 is the longitude of the 
ascending node, which defines the angle of the axis of inclination 
on the sky. 

We applied this method to all of the ellipse fits within 
la of the best-fit values listed in Table 2. Table 3 lists the 
resulting deprojected values. We constrain the disk inclination 
to 28? 5 zb 2°, consistent with previous estimates (Schneider 



Figure 5. Best-fit ellipse to the belt maximum. The projected major axis and 
true major axis are shown in red and yellow, respectively, with periastron at 
point “p.” The stellar location is marked with a yellow star. 

(A color version of this figure is available in the online journal.) 

et al. 2006), and marginally constrain the eccentricity to a non- 
zero value of 0.02 d= 0.01. This small true eccentricity results 
in a poorly constrained co. As sl result, the pericenter of the 
HD 181327 disk could be located anywhere in the SW quadrant 
of Figure 4. 

The red ellipse in Figure 5 shows the best fit to the HD 181 327 
illumination-corrected belt maximum. The red line illustrates 
the major axis of the projected ellipse and the yellow line marks 
the true major axis with periastron marked with a “p The 
line of nodes (the axis of inclination of the disk) is nearly 
coincident with the projected major axis, as expected for a nearly 
circular ring; the line connecting forward to back scattering is 
approximately perpendicular to the projected major axis. The 
center of the ellipse and the star are marked with a red dot and 
a yellow star, respectively. 

Table 3 also lists the deprojected ellipse parameters for fits to 
the inner edge of the belt. In the case of an infinitely thin and 
uniform disk, we expect ellipses fit to the belt maximum and the 
belt inner edge to agree within their mutual uncertainties, except 
for a. A discrepancy would suggest that the radial distribution 
of dust varies significantly, that the disk is not flat, or that the 
disk has an opening angle of more than a few degrees. As shown 
in Table 3, fits to the illumination-corrected belt maximum and 
the illumination-corrected inner edge give consistent values for 
all parameters and have similar uncertainties; a flat, thin disk 
appears to be a reasonable approximation for the HD 181327 
debris belt. 
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Figure 6. Deprojection procedure. Panel (a): original STIS image. Panel (b): (1/r 2 ) map normalized to the birth ring distance of 90.5 AU. Panel (c): panel (a) divided 
by panel (b). Panel (d): map of the average scattering angle. Panel (e): map of the average semi-major axis. Panel (f): best-fit scattering phase function map. Panel (g): 
panel (c) divided by panel (f). Panel (h): deprojected optical depth with the belt’s pericenter located directly to the right of the star. The stellar location is marked in 
each image with a white star. 

(A color version of this figure is available in the online journal.) 


3.3. Illumination Correction 

We next corrected for the 1/r 2 stellar illumination factor. We 
distributed 10 7 particles uniformly in azimuth and logarithmi- 
cally in r from 50 to 800 AU. Using the deprojected ellipse 
parameters from the first line of Table 3, we then rotated the 
three-dimensional positions by co , i, and Q, and binned the par- 
ticles into pixels based on their on-sky (. x , y ) coordinates, the 
pixel scale of 50.77 mas, and a distance to HD 1 81327 of 5 1 .8 pc. 
Finally, we calculated the average 1/r 2 value in each pixel to 
produce a (1/r 2 ) map, which we divided into the HD 181327 
STIS image. 

Figure 6 shows this portion of the deprojection process. Panel 
(a) shows the initial HD 181327 STIS image. Panel (b) shows 
the (1/r 2 ) map for the deprojected ellipse fit listed in the first 
line of Table 3, normalized to the radial peak at 90.5 AU. 
Panel (c) shows the illumination-corrected image, given by the 
HD 181327 STIS image divided by the (1/r 2 ) map. 


g values much less than is expected for micron sized grains 
(e.g., Kalas et al. 2005; Schneider et al. 2006; Debes et al. 2008; 
Thalmann et al. 2011). Additionally, the SPF of the zodiacal 
dust cloud is significantly flatter near a scattering phase angle 
of 90° than predicted by a single forward- scattering HG phase 
function (Hong 1985). 

Instead of using an analytic phase function, we fit an empirical 
SPF to the data. At a given semi-major axis, a , the illumination- 
corrected flux, F', shown in Figure 6(c) is proportional to 

, f dN(a, 0) 9 

F(a,6) oc / — — s 2 2sca0) pio, s)ds, (2) 

where dN(a, 0)/ds is the differential number of grains of size s, 
GscaW is the scattering efficiency, and p(0, s ) is the SPF. In the 
case of a uniform, unperturbed disk with small eccentricity, the 
size distribution is constant along a given semi-major axis, i.e., 
independent of 0. As a result, the SPF averaged over grain size 


3.4. The Scattering Phase Functions ) 
and Density Distribution 

The SPFs of disks are commonly approximated using a 
Henyey-Greenstein (HG) SPF, given by 

1 1 — g 2 

P ^ An [ 1 + g 2 — 2 g cos 0 ] 3 / 2 ’ 

where 0 is the scattering phase angle and g is the HG asymmetry 
parameter, ranging from — 1 for perfect back scattering to 1 for 
perfect forward scattering. Small dust grains are known to be 
forward scattering, so typically g > 0 for debris disks. However, 
this function is typically used out of expedience, not fidelity. 
SPFs predicted by Mie theory do not resemble HG functions in 
many cases, and HG SPF fits to observed debris disks produce 


p(a,0) 


f s 2 2sca(s) p(0, s)ds 

FWI 1 QMds 


( 3 ) 


is a function of a and 0 . Therefore we can determine an empir- 
ical SPF by fitting the illumination-corrected flux as a function 
of a and 9 . 

Following the same procedure described above for the (1/r 2 ) 
map, we created projected maps of the average scattering phase 
angle (0) (Figure 6(d)) and the average semi-major axis {a) 
(Figure 6(e)). For a given value of a, we then selected those 
pixels with \(a) — a\ < 2.63 AU (1 pixel width) and fit the 
illumination-corrected flux as a function of {0). 

Figure 7 shows the illumination-corrected flux as a function 
of (0) for two values of a. Here we used the line joining the 
minimum and maximum scattering phase angles in Figure 6 
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Figure 7. Illumination-corrected flux as a function of scattering phase angle 
at and exterior to the belt maximum (lower and upper panels, respectively). 
The SE and NW halves of the disk are shown in red and black, respectively. 
The yellow line shows the best fourth degree polynomial fit to the scattering 
phase function. The dashed green line shows the best-fit Henyey-Greenstein 
phase function, a poor fit to the observed variation. 

(A color version of this figure is available in the online journal.) 

to divide the disk into a SE half (shown in red) and a NW 
half (shown in black). In the case of a uniform disk, the SE 
and NW SPFs should be identical and this is approximately 
true at a = 105 AU, as shown in the top panel of Figure 7. 
However, the bottom panel of Figure 7 shows that near the belt 
maximum, a = 89.4 AU, the two halves are asymmetric. Under 
the assumption of a flat, thin disk, such asymmetries cannot be 
explained by additional projection or scattering effects and form 
the foundation by which we identify density asymmetries. 

We fit the flux measurements as a function of scattering phase 
angle with a fourth degree polynomial, shown in yellow in 
Figure 7, to obtain an empirical SPF. This function could be 
fit to the flux from the SE half of the disk, the NW half of the 
disk, or both halves simultaneously. In the case of HD 181327 
the SE disk flux is consistently less than the NW flux. We chose 
to interpret any detected asymmetries as density enhancements, 
so we fit the empirical SPF to the flux from the SE half of the 
disk only. 

Repeating this procedure for a ranging from 7 1 AU to 263 AU 
in steps of A a = 2.63 AU (1 pixel), we obtained the SPFs 
shown in Figure 8. Because grains beyond the parent body 
ring should be size-sorted (see Section 4.1.2), with s decreasing 
with increasing a , the normalization of the SPF at a given a is 
degenerate with the average Q sca and surface density at a given a. 
Additionally, the normalization of a given SPF strongly depends 
on the behavior of the SPF at small 0 , but our observations are 
limited to 60° < 6 < 120° given the disk inclination ~30°. 
Thus, we chose to normalize each SPF at 6 = 90°, which 
maintains the radial profile along the 0 = 90° line. We note 
that, in the end, we present any detected density asymmetries in 



Scattering angle (°) 

Figure 8. Empirically derived scattering phase function as a function of a. The 
functional behavior is consistent with Mie theory predictions for smaller grains 
at larger circumstellar distances. 

(A color version of this figure is available in the online journal.) 


terms of their detection confidence, which is independent of the 
normalization. 

Using the empirical SPFs and the (0) maps shown in 
Figure 6(d), we created the SPF map shown in Figure 6(f). We 
divided the illumination-corrected surface brightness image by 
the scattering phase function map to produce the projected op- 
tical depth shown in Figure 6(g). Finally, rotating by Q to align 
the longitude of nodes along the x axis, stretching the image 
vertically by 1 / cos i , and rotating by co to place the periastron 
on the positive x axis gives the face-on optical depth shown in 
Figure 6(h). 

Figure 9 shows “radial” profiles for the HD 181327 disk, 
over the region for which we are confident PSF residuals 
are insignificant. The left panel shows cuts along the major 
and minor axes of the illumination-corrected flux density as a 
function of semi-major axis using a median binning with width 
Aa = 1 pixel. The FWHM of the illumination-corrected flux 
varies depending on the cut, from ~25% to ~40%. Normalizing 
the SPF to a scattering angle of 90° produces an optical depth 
profile with FWHM of 30%, while normalizing to a scattering 
angle of 60° produces an optical depth profile with FWHM of 
40%. We do not know the correct normalization for the empirical 
SPF as a function of a, so the true radial dependence and FWHM 
of the optical depth is unknown. We attempted to normalize the 
SPFs by extrapolation to small scattering angles using a variety 
of appropriate functions; the resulting radial power law of the 
optical depth was wildly uncertain. Thus, while the FWHM may 
help constrain the mass of a planet sculpting the inner edge of 
the disk in the case of a single SPF (Rodigas et al. 2014), if the 
SPF varies significantly with radial distance these constraints 
are less certain. 

The right panel shows the median radial profile of the 
deprojected normalized optical depth when normalizing the SPF 
to a scattering angle of 90°. The inner edge of the optical 
depth is oca 7,0 . The power law of the outer edge varies 
smoothly with semi-major axis. From 95-140 AU, the optical 
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Figure 9. Left: semi-major axis profile of the illumination-corrected flux density, median binned using a width Aa = 1 pixel, along the apparent major and minor 
axes. Right: median semi-major axis profile of the deprojected optical depth shown in panel (h) of Figure 6 (solid black line) along with three power-law fits (colored), 
where ao = 90.5 AU. The dotted line shows the deprojected optical depth profile when normalizing the empirical SPF to a scattering angle of 60°; the FWHM of the 
radial optical depth profile is degenerate with a radial-varying SPF. 

(A color version of this figure is available in the online journal.) 
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Figure 10. Optical depth deviations from a uniform disk. Left: the on-sky projected plane. Middle: smoothed residuals in the deprojected face-on plane. Right: 
detection confidence of the smoothed, deprojected residuals. 

(A color version of this figure is available in the online journal.) 


depth aa -3 ' 7 , while from 150-250 AU the optical depth is 
oca~ lJ . The dotted line shows the median radial profile of the 
deprojected normalized optical depth when normalizing the SPF 
to a scattering angle of 60°. The radial optical depth profile is 
degenerate with a radially dependent SPF. 

3.5. Detected Asymmetries 

The normalized optical depth in Figure 6(h) shows one 
obvious asymmetry: a density enhancement in the birth ring 
extending ~90° counter-clockwise from periastron. To reveal 
other less evident asymmetries, we took the difference between 
Figure 6(h) and its smooth counterpart, i.e., a uniform eccentric 
disk. To create the deprojected optical depth for a uniform 


eccentric disk, we could simply calculate the median of the 
deprojected optical depth as a function of a. However, we chose 
to interpret any asymmetries as density enhancements, so we fit 
the optical depth as a function of a using a smoothly varying 
polynomial, then set the optical depth of the smooth disk equal 
to the minimum of this function. 

Figure 10 shows the fractional optical depth residuals from 
what would be expected for a flat, thin, uniform disk. The left 
panel shows these normalized residuals in the projected sky 
plane. The middle panel shows the same residuals smoothed 
using a 3 x 3 median boxcar in the deprojected, face-on plane 
with periastron located along the positive v axis. The right panel 
shows the detection confidence of the smoothed, deprojected 
residuals. 
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4. DISCUSSION 

We have detected asymmetric residuals in the HD 181327 
debris disk at >3 o over many AU, suggesting the possibility of 
density enhancements/deficits. However, to derive the residuals 
shown in Figure 10, we assumed the disk was thin and flat. Thus, 
there are three possible causes for the detected asymmetries. 
First, the disk could indeed be thin and flat, implying that 
the detected asymmetries reflect true density enhancements or 
deficits. Second, our assumption that the HD 181327 disk is thin 
could be incorrect. Third, our assumption that the HD 181327 
disk is flat could be incorrect. We discuss each of these scenarios. 

4.1. Density Enhancements in a Thin, Flat Disk 

In the case of a thin, flat disk, the observed asymmetries can- 
not be explained by projection or scattering effects; the detected 
asymmetries necessitate density asymmetries. We note that we 
chose to interpret the asymmetries as density enhancements. An 
equally valid, but less likely explanation is that the asymmetries 
represent deficits of material elsewhere in the disk. Similarly, 
the detection confidence in the right panel of Figure 10 rep- 
resents the confidence of an asymmetry, not the confidence of 
whether an asymmetry is an excess as we have assumed, or 
deficit elsewhere in the disk. 

4.1.1. A Massive Collisional Event 

Figure 10 shows that the density enhancement appears as 
an arc of material near periastron in the birth ring extending 
counter-clockwise and radially outward, with its angular size 
increasing with circumstellar distance. The geometry of the 
detected asymmetry is consistent with models of small grains 
produced by high-mass collisional events (e.g., Grigorieva et al. 
2007; Jackson & Wyatt 2012; Krai et al. 2013; Jackson et al. 
2014). In these models, the smallest fragments with > 0.5 
quickly leave the system on hyperbolic orbits. However, the 
largest fragments, unaffected by radiation pressure, remain 
bound to the system on orbits that all cross at the location of 
the initial impact, creating a natural “pinch” point. These large 
grains routinely collide at this “pinch,” producing a continual 
outflow of small grains. The lifetime of this asymmetry varies 
dramatically among models, from a few orbital periods to 
millions of years, and is largely dependent on the fraction of 
the dust mass associated with the presumed massive collision. 

Interpreting the detected asymmetry as originating from a 
collisional event, we can place a lower limit on the mass of dust 
generated by the collision. We multiplied the middle panel of 
Figure 10 by the optical depth peak value of ro = 2 x 10 -3 , 
estimated from Li R /L* for HD 181327 (Lebreton et al. 2012), 
to obtain the absolute size-integrated optical depth. We then 
multiplied by the pixel area, calculated from the observed 
image scale of 50.77 mas pixel -1 , and an assumed distance 
to HD 181327 of 51.8 pc (Holmberg et al. 2009) to obtain 
the absolute size-integrated product of the cross section and 
scattering efficiency, g sca . We summed these values for all 
pixels with an S/N > 2cr to calculate a total <2 sca -weighted 
cross section associated with the collisional event. We assumed 
the particles are well-described as porous mixtures of ice, 
amorphous silicate, and carbonaceous material, as determined 
from Mie theory by Lebreton et al. (2012), and calculated 
the mass associated with the total <2 sca -weighted cross section, 
where g sca was averaged over the STIS bandpass. We assumed a 
size distribution dN /ds oc a - 3 85 , the expected size distribution 
of fragments from a recent disruptive collision (Leinhardt & 


Stewart 2012), though the results are relatively insensitive to the 
size distribution assumed because Q SC3L peaks near s = l /xm. 

We find that the detected asymmetry requires 10 20 kg of dust 
smaller than a few microns, equivalent to 1 % the mass of Pluto or 
~10 -4 the mass of the Kuiper Belt. However, we only examined 
the asymmetric component of the optical depth and may have 
removed a significant fraction of the dust associated with this 
collisional event. Therefore, the estimated dust mass represents 
a lower limit. 

Calculating the total mass involved in the collision requires 
extrapolating the estimated mass over a narrow range of micron- 
sized dust grains to bodies kilometers in size. The power law 
used for such an extrapolation is not well understood and may 
feature a number of breaks. We therefore examine the lower 
limit on the total collisional mass, i.e., 100% of the target mass 
is converted into grains smaller than a few microns. Given the 
target mass of 10 20 kg, the binding energy of the target can be 
approximated as the gravitational binding energy, or 10 24 J. The 
collisional energy is given by 

£coi = ft;, 2 , (4) 

where /x = m t m p /(m t + m p ), m p is the projectile mass, m t is 
the target mass, and v t is the impact velocity. Given the disk’s 
small eccentricity, we assume the impact velocity is dominated 
by the scale height of the disk and set v t = (H/r) i; or bit, where 
t? orbit = 3.5 km s -1 is the orbital speed in the parent ring, 
and EL /r = 0.1, the maximum value allowed by our analysis 
(see Section 4.2). Assuming half of the collisional energy is 
imparted to the target and 100% of the imparted energy goes 
toward fragmenting the target, we estimate ra p ~ 10 20 kg, i.e., 
the projectile and target mass are roughly the same mass. 

Explaining the observed dust excess using a single collisional 
event becomes difficult when increasing the target mass beyond 
the lower limit of 10 20 kg. For example, assuming a size distribu- 
tion dN /ds oc a -3 - 85 valid for all sizes, we estimate a total colli- 
sional mass of 10 22 kg, roughly the mass of Pluto. To catastroph- 
ically fragment a Pluto mass object, we must increase the impact 
velocity by 750 m s -1 . We have ruled out values of H/r > 0.11, 
so we must invoke more exotic scenarios, like a massive com- 
pact binary target or an unseen planet stirring the disk, to explain 
the catastrophic disruption of a Pluto mass object. 

Alternatively, one could argue that the observed dust excess 
is not the result of a single recent collision, but a collisional 
avalanche initiated by a much less massive target. Grigorieva 
et al. (2007) showed that the initial release of 10 17 kg of dust in 
a disk with r ~ 10 -3 can trigger a collisional avalanche with a 
total dust cross section that peaks at 200 times the initial cross 
section released. Models of collisional avalanches (Grigorieva 
et al. 2007; Krai et al. 2013) also predict morphologies qualita- 
tively consistent with the observations. A more detailed model 
of a collisional avalanche in the HD 181327 disk is required to 
determine the probability of witnessing such an event. 

4.1.2. The SPF and Size Segregation 

An unperturbed narrow ring of parent bodies producing dust 
should naturally produce a size- sorted halo of dust exterior to the 
parent body ring. Upon launch, a dust grain’s orbit is modified 
by radiation and solar wind pressure. The post-launch apastron 
distance of a dust grain increases with /3, where /3 is the ratio 
of radiation force to gravitational force on the dust grain. For 
many materials and stellar systems, /3 oc s~ l is approximately 
valid all the way down to the blowout grain size. Smaller grains 
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Figure 11. SPF predicted by Mie theory for the best-fit composition of Lebreton 
et al. (2012). The trend suggests grain size decreases with increasing semi-major 
axis in the HD 181 327 debris disk, as predicted by dynamical debris disk models. 
(A color version of this figure is available in the online journal.) 



Scattering angle (°) 

Figure 12. Best- fit SPF predicted by Mie theory for the best-fit composition of 
Lebreton et al. (2012). These fits require sub-blowout-size grains to dominate 
the HD 181327 optical depth. 

(A color version of this figure is available in the online journal.) 


achieve larger apastron distances and dominate the cross section 
density at larger circumstellar distances. 

The empirical SPFs in Figure 8 are consistent with this size 
sorting (orange through purple curves). These curves show that 
over the range of observed scattering angles, the width of the 
forward-scattering peak increases with semi-major axis. This 
trend is consistent with what is expected from Mie theory if the 
dominant grain size decreases with increasing semi-major axis. 

To illustrate how the empirical SPF suggests size segregation, 
Figure 1 1 shows the normalized SPF predicted by Mie theory as 
a function of grain size for the porous mixture of ice, amorphous 
silicate, and carbonaceous material determined from fits to the 
spectral energy distribution (SED) of HD 181327 as suggested 
by Lebreton et al. (2012). Although the magnitude of forward 
scattering for these grains does not agree with that shown in 
Figure 8, the trend suggests size segregation in the HD 181327 
debris disk with particle size decreasing as semi-major axis 
increases. This trend is a broad prediction of Mie theory; other 
compositions show similar trends. 

While models show that larger grains should dominate the 
halo’s optical depth closer to the birth ring, these same models 
predict that the trend stops at the outer edge of the birth ring. 
As shown by Strubbe & Chiang (2006) and Thebault et al. 
(2014), the dominant grain size within the birth ring should 
be approximately equal to the blowout size, i.e., the smallest 
grains in the system. If this were the case for HD 181327, we 
would expect to see the trend in the empirical SPF reverse upon 
reaching the birth ring, i.e., the SPF at a = 89.4 AU should be 
similar to that at a = 210 AU. Instead, the empirical SPF in the 
parent body ring (black and red curves in Figure 8) continues 
the observed trend, with an apparent decrease in the degree of 
forward scattering of the range of observed scattering angles. 
Additionally, the degree of back scattering appears to increase 
and the SPF flattens near a scattering phase angle of 90°. 

If the empirical SPF in the birth ring closely matches the true 
SPF, then there may be an absence of small grains in the birth 
ring. Such a scenario is expected if a planetary companion orbits 
exterior to the birth ring. As shown by Thebault et al. (2014), 


gravitational perturbations from an exterior planet can dynami- 
cally eject small grains, whose orbits are planet-crossing, before 
they can appreciably contribute to the optical depth at their peri- 
astron distance (the birth ring), while leaving the expected size 
sorting trend beyond the orbit of the planet unaffected. 

The degree of forward scattering predicted by Mie theory, as 
shown in Figure 1 1 , does not closely match the empirical SPF 
in the birth ring, even for very large grains. Is the empirical 
SPF measured in the birth ring reasonable for large debris disk 
grains? The thin dashed line in Figure 8 shows the zodiacal 
cloud’s derived SPF (assuming v = 1, see Hong 1985), which 
is dominated by ~100 /xm grains near 1 AU (Grun et al. 
1985). Over the range of observed scattering angles, the zodiacal 
cloud’s SPF is less forward scattering than that observed for 
HD 181327; the SPF from 100 /xm zodiacal cloud grains is also 
inconsistent with Mie theory, but consistent with the low degree 
of forward scattering observed in HD 181327. 

Alternatively, the SPF predicted by Mie theory for ~0.1 /xm 
grains fits the empirical SPF in the birth ring, and grains that 
increase in size with semi-major axis can fit the empirical SPF all 
the way out to 210 AU, as shown in Figure 12. However, such 
small grains are well below the blowout size for this system 
and do not survive long enough to dominate over the bound 
grains. Further, for an unperturbed debris disk there is no known 
physical mechanism to create the trend of increasing grain size 
with semi-major axis shown in Figure 12. Interestingly, the best- 
fit SED model of Lebreton et al. (2012) also requires grains 
below the blowout size, ~0.8 /xm. There are three ways to 
interpret these results. First, the best-fit porous composition 
of Lebreton et al. (2012) is correct, the optical properties of 
the HD 181327 dust are truly dominated by grains below the 
blowout size, and we have a poor understanding of the dynamical 
behavior of such grains. Second, the SPF predicted by Mie 
theory is not valid for the complex porous grains of Lebreton 
et al. (2012). Third, the SED modeling of Lebreton et al. (2012) 
needs to be revisited. Given that Lebreton et al. (2012) assumed 
a single power-law size distribution valid at all points in the 
disk, we suspect the truth is a combination of all of these points. 
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Figure 13. Extrapolated fits to the SE illumination-corrected flux as a function 
of scattering phase angle at the location of the belt maximum. The data is 
shown in red (SE) and black (NW). The solid orange line shows the best- fit 
two-component HG SPF and the dashed green line shows best-fit single HG 
with g = 0.3. The apparently low albedo of the HD 181327 disk may be entirely 
explained by a more strongly forward-scattering phase function (orange curve), 
though the extrapolation is not well- constrained. 

(A color version of this figure is available in the online journal.) 

All of the empirical SPFs shown in Figure 8 deviate signif- 
icantly from HG functions. The dashed green lines in Figure 7 
show the best-fit HG SPFs for a = 89.4 and a = 105 AU. In 
each case, the slope of the empirical SPF is significantly larger 
than the best-fit HG SPF at 0 ~ 60°. Extrapolating this slope 
to even smaller, unobservable scattering angles could suggest 
that the HD 181327 disk is far more forward scattering than 
previously thought. 

The enhanced forward scattering that we find may help 
explain the apparently low albedo for the HD 181327 disk. 
Assuming a single HG SPF with g = 0.3, Lebreton et al. (2012) 
showed that the observed albedo for the HD 181 327 disk is lower 
than their model predictions by a factor of ~4. We extrapolated 
the empirical SPF in the birth ring to all values of 0 using a 
two-component HG fit. We found that two HG SPFs with 
g% = 0.87 and g 2 = —0.30, weighted at 87% and 13%, 
respectively, best fit the illumination-corrected SE flux (red 
curve in the lower panel of Figure 7), with a x 2 / v = 2.2. The 
resulting SPF, shown in Figure 13, is significantly more forward 
scattering than a HG SPF with g = 0.3. As a consequence, the 
normalized empirical SPF is approximately 25% that of a HG 
SPF with g = 0.3 at 0 = 90°, eliminating the discrepancy noted 
in Lebreton et al. (2012). 

Unfortunately the limited range of observable scattering 
phase angles prevents us from significantly constraining the 
empirical phase function at small scattering angles. Considering 
all two-component HG fits to the SE flux with x 2 / v within a 
factor of two of the best- fit acceptable, values of g\ = [0.3 , 1.0]; 
the extrapolated SPF is statistically consistent with the albedo 
discrepancy noted in Lebreton et al. (2012). However, if the 
observed increase in flux at small scattering angles is real and 
due to the true SPF, we regard fits with low g \ with skepticism, 
as they do a poor job reproducing the SE flux at small scattering 
angles. 

The discrepancy between the observed and modeled albedo 
from Lebreton et al. (2012) may also be an artifact of size seg- 
regation. As previously noted, Lebreton et al. (2012) modeled 
the disk using a single grain size distribution at all circumstellar 
distances. Instead, let us suppose the size distribution resembles 


that implied by the empirical SPF, with a paucity of small grains 
in the birth ring, and a halo dominated by small grains that de- 
crease in size at larger circumstellar distances. In this case the 
most efficient scatterers, the sub-micron grains, would not con- 
tribute as significantly to the scattered light image. The thermal 
emission, however, will be dominated by the large grains in the 
parent body belt closer to the star. As a result, the observed 
albedo will be lower than in the uniform size distribution case. 

4.1.3. Alternative Density Distributions 

We fit the illumination-corrected surface brightness 
of the SE half of the disk with an empirical SPF and discussed 
the implications of such an SPF above. However, one could ar- 
gue that the empirical SPF does not reflect the true SPF, and the 
SPF is degenerate with some non-uniform density distribution. 
While strictly true, any observed surface brightness distribu- 
tion can be attributed to a contrived density distribution. For 
example, additional density enhancements near scattering an- 
gles of 90° could cause the flatness of the empirical SPF in the 
birth ring. In light of this, we limit ourselves to two relevant 
physical scenarios: (1) a birth ring that is in fact occupied by 
small dust grains but has a density enhancement that masks the 
true SPF, and (2) a debris disk whose true SPF everywhere is 
approximately equal to the empirical SPF in the birth ring. 

To investigate the nature of the necessary density enhance- 
ments we deprojected the disk again, but altered the SPF in 
panel (f) of Figure 6. For the first case, in which we assume 
the birth ring is occupied by small dust grains, we assigned the 
empirical SPF from pixels with 208 < a < 213 AU to all pixels 
with 86.8 < a < 100. For the second case, we assigned the 
empirical SPF from pixels with 86.8 < a < 92.0 AU to all 
other pixels. 

Figure 14 shows the results of both of these scenarios in the 
deprojected frame with periastron pointing to the right. For the 
true SPF in the birth ring to match the empirical SPF in the outer 
disk, as expected for an unperturbed narrow birth ring, there 
must be an additional density enhancement in the birth ring. 
The left panels show that this density enhancement is ~100% 
near periastron and is superimposed on top of the asymmetry 
previously detected. The additional density enhancement must 
be symmetric and aligned, by chance, perpendicular to the line 
of nodes, a scenario we find unlikely. 

If we instead assign the empirical SPF of the birth ring to all 
other points in the disk, the required density enhancement in the 
outer disk increases significantly, as shown in the right panels 
of Figure 14. Essentially the density distribution must make up 
for the lack of forward scattering in the birth ring’s SPF. This 
interpretation requires a slightly more massive collisional event, 
with 7.5 x 10 20 kg of dust less than a few microns in size, or the 
equivalent of 7.5% the mass of Pluto. 

4.2. Constraints on the Scale Height 
of the HD 181327 Debris Disk 

A non-zero scale height, which we have so far ignored, 
would affect the appearance of a debris disk in two ways. 
First, for a disk with a vertical density distribution peaked at 
the midplane, the line of sight intersects more dust near the 
midplane when looking near the ansae; a non-zero scale height 
would brighten the disk along the ansae. Second, a non-zero 
scale height would change the amount of flux received along 
the line joining forward to back scattering. Along this axis, the 
line of sight intercepts a wider range of circumstellar distances 
and scattering angles. 
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Small Grains in Birth Ring Birth Ring SPF Applied Everywhere 



Figure 14. Alternative deprojected optical depth distributions. Left: the required optical depth distribution if small grains dominate both the birth ring and outer disk, 
as expected from unperturbed disk theory (obtained by setting the SPF of the birth ring equal to the SPF of the outer disk). Right: the required optical depth distribution 
if the empirical SPF in the birth ring is true everywhere. Both scenarios feature a greater degree of asymmetry than our preferred deprojection shown in Figure 10. 

(A color version of this figure is available in the online journal.) 


To constrain the scale height of HD 181327, we examined a 
sharp radial feature that a non-zero scale height would tend to 
blur: the inner edge of the observed belt. We fit a radial power 
law to the belt’s flux as a function of projected semi-major axis 
at four locations within the inner edge (71.0 < a < 86.8 AU): 
the east and west disk ansae (along the axis of inclination), 
along the line of forward scattering, and along the line of back 
scattering. The power law at the ansae should agree in the case of 
a uniform disk. Unfortunately our observations are least certain 
in these regions, so we averaged the two power laws together. We 
found radial power-law exponents of 4.5 zb 0.4, 4.3 =b 0.3, and 
4.0 d= 0.5 for the ansae, forward- scattering, and back-scattering 
sides of the disk, respectively. 

We then produced Monte Carlo models of circularly sym- 
metric disks including single- scattering radiative transfer for 
comparison. We assumed a “knife-edge” radial density profile, 
with n(r) a r & for r < r pea k, and n(r) a r a for r > r pea k- We 
investigated values of a in the range [—12, —2] and in the 
range [4, 10], and 88.1 ^ r pea k ^ 94.7 AU. We investigated 
disk scale heights 0 < H/r <0.3 and used a HG SPF to rep- 
resent the SPF of the disk with 0 < g < 0.9. We produced a 
total of 1.6 million models, each using 1 million particles, and 
convolved each model’s image with a Tiny Tim PSF, based on 
the stellar spectral template for an F6, B — V = +0.42 star. 


We calculated the power laws at the inner edges of our models. 
We found that only models with H/r < 0.11 produced profiles 
that were simultaneously within 1 o of the measured values at the 
ansae, forward-scattering, and back- scattering sides of the disk. 
Our constraint is consistent with the H/r < 0.09 constraint 
from Schneider et al. (2006). 

Given our constraint on the scale height of the disk, could 
a non-zero scale height explain the variation in the empirical 
SPF observed for HD 181327? To check, we calculated the 
apparent SPF for each of our Monte Carlo Models. We found that 
only scale heights >0.3 produced a significant decrease in the 
apparent forward scattering of the birth ring SPF, inconsistent 
with the constraints on the scale height. Further, we were unable 
to qualitatively reproduce the curves shown in Figure 8 from 
our models, even for large scale heights. Figure 15 shows an 
example of our model results for H/ r = 0.3. The models cannot 
produce SPFs with the observed spread in forward scattering 
from circumstellar distances of 105-210 AU. Most critically, the 
models produced SPFs that were nearly tangential at 0 = 90°. 
This is because a thick disk with a single radial power law 
beyond the birth ring reduces both forward and back scattering 
close to the disk, a trend that is not observed in HD 181327. We 
conclude that a non-zero scale height alone cannot explain the 
variation in the empirical SPF. 
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Figure 15. SPF vs. a for our Monte Carlo model with a scale height of 
H/r = 0.3. A non-zero scale height decreases forward scattering in the birth 
ring as observed, but cannot produce the shape or spread of the empirical SPF. 
(A color version of this figure is available in the online journal.) 

The radial flux distribution in the outer disk is also inconsis- 
tent with the scenario of a thick disk with a single forward SPF. 
Power-law fits to the flux,/, in the region 105 < a < 131 AU 
reveal / a a~ 5 6 along the ansae, but / a a -51 in the direction 
of forward scattering and / a a~ 6 9 in the direction of back 
scattering. This trend is found regardless of the inner and outer 
boundaries of the fit as long as a > 90.5 AU, and is robust 
within the range of acceptable ellipse inclinations. Our Monte 
Carlo models reveal that a thick disk with a single forward SPF 
should have a radial power law that is steeper in the direction of 
forward scattering, not shallower. 

Conversely, the shallower profile in the direction of forward 
scattering is consistent with a thin disk with an ^-dependent 
SPF. To illustrate this, we used our Monte Carlo code to 
produce a simple circularly symmetric model of a thin disk 
with cross section density ocr -2 . We assigned each particle in 
the model a HG SPF with g = 0.9(r/158 AU). In the region of 
105 < r < 131 AU, we found / oc r 5 6 along the ansae while 
/ a r -5 0 along the line of forward scattering. 

We conclude that the HD 181327 debris disk scale height is 
not sufficiently large to produce any of the observed asymme- 
tries or trends in the SPF. We therefore consider the HD 181327 
disk to be “thin.” 

4.3. Possible Warping of the HD 181327 
Debris Disk: ISM Interactions 

A number of debris disks exhibit clear signatures of ISM inter- 
actions (e.g., Hines et al. 2007; Debes et al. 2009). Interactions 
with the ISM gas can force dust grains out of the disk midplane, 
warping an otherwise thin, flat disk into three dimensions. Our 
deprojection assumed the disk was thin and flat, such that the 
true scattering angle was constant along a given radial line to the 
star. Warping by ISM gas could effectively cause the scattering 
angle to change with circumstellar distance, creating projection 
and scattering effects that were not previously considered. 

To test whether ISM gas could cause the asymmetries detected 
in the STIS observations of the HD 181327 debris disk, as 
well as the observed trend in the empirical SPF, we produced 
simple models of debris disks interacting with ISM gas. We 


considered only grains below the blowout size. Assuming the 
best-fit dust grain composition determined by Lebreton et al. 
(2012), moderate ISM interactions should result in grains below 
the blowout size dominating the STIS observations of the 
HD 181327 debris disk. We calculate that in the absence of 
ISM interactions, blowout grains contribute ~ 10% of the optical 
depth in the parent ring over the STIS bandpass, in spite of their 
short life times. Barely-bound grains, predicted to dominate the 
parent ring cross section in the absence of ISM interactions 
(Strubbe & Chiang 2006) should be easily entrained by the ISM 
at apastron and removed from the system quickly, resulting in 
the dominance of blowout grains. 

We launched 200,000 dust grains from a 10 AU wide uniform 
parent ring centered at 90 AU. We modeled 30 dust grain sizes 
below the blowout size, spaced logarithmically from 0.1 to 
4.9 pi m. We used the best-fit dust grain composition obtained 
by Lebreton et al. (2012) and converted grain size to /3 using 
Figure 12 in Lebreton et al. (2012). We integrated the equations 
of motion described by Debes et al. (2009) using a fourth order 
Runge-Kutta method and a time step size equal to 0.05% 
the orbital period at the birth ring. We integrated until the 
stellocentric distances of all grains exceeded 10 times their 
initial stellocentric distances. We considered both prograde and 
retrograde orientations for the orbits. 

HD 181327 has a proper motion of 23.99 mas yr -1 in right 
ascension and —81.82 mas yr -1 in declination (van Leeuwen 
2007) and has a negligible radial velocity (Gontcharov 2006), 
implying a velocity vector in the plane of the sky of 21 km s -1 
oriented at 16? 3 E of S given a distance of 51.8 pc. We 
examined 154 relative velocity vectors v re i = r* — vism between 
HD 181327 and the ISM gas, covering seven values of speed 
from 1 to 100 km s -1 and 22 directions. We set the density of 
the ISM gas to 1.67 x 10 -22 g cm -3 , the value used by Debes 
et al. (2009), and used a stellar mass of 1 .36 M© (Lebreton et al. 
2012). 

We used dustmap to produce images of the disk (Stark 2011) 
using the size distribution and gsca values obtained by Mie the- 
ory for the best-fit model of Lebreton et al. (2012). We assumed 
a single HG SPF valid for all grain sizes. We examined seven 
values of the SPF asymmetry parameter g ranging from 0.2 to 
0.9. We note that while the observed SPF for HD 181327 does 
not resemble a HG SPF, we chose to use the simple function for 
numerical rapidity; we do not intend to find a quantitative best fit 
to the observations. We then reduced the modeled data using the 
same methods described above and searched for asymmetries 
qualitatively similar to those detected in the STIS image. 

The top row of Figure 16 shows an example model in which 
t> rel = 50 km s -1 . We show the STIS observations of HD 181327 
in the bottom row for direct comparison. The model residuals 
(Figure 16(c)) appear qualitatively similar to those observed in 
the STIS image (Figure 16(f)). We are able to produce similar 
structures for both prograde and retrograde orbits. We find that 
ISM interactions can qualitatively reproduce the bright arc along 
the S W portion of the parent ring if the relative velocity between 
HD 181327 and the ISM exceeds ~30 km s -1 and is oriented 
toward the SW. 

Given a minimum relative velocity iy e] ~ 30 km s -1 oriented 
toward the SW, and a stellar velocity of 21 km s -1 at 16? 3 E of 
S for HD 181327, the ISM velocity must contribute strongly to 
the total relative velocity. Roughly speaking, viSM must have 
an eastward component >25 km s _1 . This is in contrast to other 
ISM-sculpted debris disks that exhibit geometries generally 
consistent with u re i being dominated by the stellar velocity. 
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Figure 16. Our best ISM-perturbed disk model (top row) compared to the STIS observations (bottom row), reduced with identical pipelines. Panel (a): model image. 
Panel (b): deprojected face-on optical depth of the model. Panel (c): smoothed residuals of the model in the deprojected face-on plane. Panel (d): HD 181327 STIS 
observations. Panel (e): deprojected face-on optical depth of STIS observations. Panel (f): smoothed residuals of the STIS observations in the deprojected face-on 
plane. The stellar location is marked in each image with a white star. 

(A color version of this figure is available in the online journal.) 


Strong ISM interactions, which should deplete barely-bound 
grains and increase the relative contribution of blowout grains, 
may help explain the sub-micron minimum grain size inferred 
from fits to the SED (Lebreton et al. 2012). Additionally, 
sub-micron blowout grains, which should dominate the STIS 
observations, may help explain the observed degree of forward 
scattering, as illustrated in Figure 12. 

However, our ISM-perturbed disk models have a number of 
caveats. First, blowout grains appear to be unable to explain 
the trends in the SPF as a function of semi-major axis shown 
in Figure 8. The model shown in Figure 16 was imaged using 
a HG SPF with g = 0.4 for all dust grains, independent of 
grain size. Figure 17 shows the empirical SPF derived using the 
techniques described in Section 3.4. The warping of the disk due 
to ISM interactions qualitatively leads to the right semi-major 
axis trend in the SPF for scattering angles >90°, but the models 
cannot reproduce the correct semi-major axis trend in the SPF 
at scattering angles <90° (cf. Figure 8). We also investigated 
more complex SPFs, including a size-dependent SPF calculated 
from Mie theory and the empirical SPF measured in the birth 
ring of the STIS observations, but were unable to reproduce the 
trends shown in Figure 8. 

Second, the majority of ISM-perturbed disks show clear signs 
of a bow shock. If v re i is indeed oriented toward the SW, then 
the HD 181327 disk should show signs of a bow shock in the 
SW quadrant of the STIS observations; we see no signs of such 
a feature out to circumstellar distances ~500 AU. 



Figure 17. Normalized SPF for an example ISM-perturbed disk model with 
Hsm = 40 km s -1 in the direction of 37° E of N and g = 0.4. Our ISM- 
perturbed disk models can qualitatively reproduce the correct trend in the SPF 
as a function of a for scattering angles >90°, but cannot produce the correct 
trend for scattering angles <90° (cf. Figure 8). 

(A color version of this figure is available in the online journal.) 


15 



The Astrophysical Journal, 789:58 (16pp), 2014 July 1 


Stark et al. 


Finally, to produce the asymmetry shown in Figure 16, we 
assumed an ISM density ~100 H cm -3 , roughly three orders 
of magnitude greater than observed in the local bubble (Frisch 
et al. 2009), where HD 181327 resides (Lallement et al. 2014). 
On the other hand, the 35 pc distant HD 61005 shows clear signs 
of what is thought to be ISM interaction (Hines et al. 2007), 
and does not appear to be near any localized enhancements in 
ISM opacity within the local bubble (Lallement et al. 2014). It 
is unclear why we observe signs of ISM interactions in disks 
within the local bubble; our dynamical understanding of dust 
undergoing ISM drag may need revision. Nonetheless, the large 
ISM density values required to produce asymmetries in this 
work are similar to those used for other nearby ISM-perturbed 
debris disks (Debes et al. 2009). 

5. CONCLUSIONS 

We have imaged the HD 181327 debris disk with STIS using 
six-roll PSF-template subtracted coronagraphy and processed 
it with a new multi-roll residual removal routine to further 
reduce quasi-static PSF residuals. The STIS observations reveal 
the HD 181327 debris ring in its entirety. The debris ring has 
a sharp inner edge and extended outer halo, consistent with 
a parent belt of planetesimals collisionally producing dust, 
and features a prominent azimuthal asymmetry. Using a new 
iterative deprojection procedure, we find the disk is inclined by 
28? 5 zb 2° from face-on. The parent ring density profile peaks 
at 90.5 AU from the star assuming a distance to HD 181327 of 
5 1 .8 pc, and is nearly circular (e = 0.02 d= 0.01) with periastron 
located in the SW quadrant. 

The empirical scattering phase function of the disk is non- 
Henyey-Greenstein and, due to a dramatic rise near the smallest 
observable scattering angles, appears significantly more forward 
scattering than previously thought, helping to explain the disk’s 
low apparent albedo. The empirical scattering phase function 
also varies with semi-major axis, with more distant stellocentric 
distances exhibiting a greater degree of forward scattering over 
the range of observed scattering angles. The scale height of 
the HD 181327 debris disk H/r < 0.11, such that line-of-sight 
effects due to the disk thickness are negligible, and is insufficient 
to explain the observed trends in the scattering phase function. 
If the true scattering phase function varies with semi-major axis 
as suggested by empirical fits, then the true radial profile of this 
disk’s optical depth, and possibly others’ measured to date, is 
highly uncertain. 

Assuming a flat disk, we deprojected the HD 181327 debris 
disk and removed the empirical scattering phase function to 
reveal the minimally asymmetric face-on optical depth. We 
find remaining asymmetries in the face-on optical depth. The 
morphology of these asymmetries can be explained either 
by a flat disk with density enhancements due to a massive 
collisional event, or qualitatively by a disk warped by strong 
ISM interactions. 

If the disk is flat and the asymmetry is due to a collisional 
event, the collisional mass must be greater than 10 20 kg, or 1% 
the mass of Pluto. The observed trends in the scattering phase 
function beyond the birth ring are consistent with the radial grain 
size sorting predicted by models. However, in contradiction to 
unperturbed debris ring models, the scattering phase function in 
the birth ring suggests large grains dominate the optical depth, 
possibly due to perturbations from an undetected planetary 
companion exterior to the debris ring. 

If the disk is warped by the ISM, our preliminary ISM- 
perturbed disk models suggest the relative velocity between 


the ISM and star must be oriented toward the SW, at nearly a 
right angle to the observed stellar velocity. This relative velocity 
requires an extremely dense ISM with an eastward velocity 
component >25 km s _1 . Our ISM-perturbed disk models cannot 
explain the observed changes in the degree of forward scattering 
as a function of semi-major axis and there are no signs of a bow 
shock, as is observed in other disks purportedly undergoing 
strong ISM interactions. 

Based on observations made with the NASA/ESA Hubble 
Space Telescope , from program #12228. Support for program 
#12228 was provided by NASA through a grant from the 
Space Telescope Science Institute, which is operated by the 
Association of Universities for Research in Astronomy, Inc., 
under NASA contract NAS 5-26555. C.C.S. acknowledges the 
support of a Carnegie Fellowship and an appointment to the 
NASA Postdoctoral Program at NASA Goddard Space Flight 
Center, administered by Oak Ridge Associated Universities 
through a contract with NASA. This work was also supported 
by NASA Astrobiology Institute grant NNA09DA81A. 
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